Calculating shear viscosity with confined non-equilibrium molecular dynamics: a case study on hematite - PAO-2 lubricant.

The behaviour of confined lubricants at the atomic scale as affected by the interactions at the surface-lubricant interface is relevant in a range of technological applications in areas such as the automotive industry. In this paper, by performing fully atomistic molecular dynamics, we investigate the regime where the viscosity starts to deviate from the bulk behaviour, a topic of great practical and scientific relevance. The simulations consist of setting up a shear flow by confining the lubricant between iron oxide surfaces. By using confined Non-Equilibrium Molecular Dynamics (NEMD) simulations at a pressure range of 0.1-1.0 GPa at 100 °C, we demonstrate that the film thickness of the fluid affects the behaviour of viscosity. We find that by increasing the number of lubricant molecules, we approach the viscosity value of the bulk fluid derived from previously published NEMD simulations for the same system. These changes in viscosity occurred at film thicknesses ranging from 10.12 to 55.93 Å. The viscosity deviations at different pressures between the system with the greatest number of lubricant molecules and the bulk simulations varied from -16% to 41%. The choice of the utilized force field for treating the atomic interactions was also investigated.

A Workflow for Identifying Viable Crystal Structures with Partially Occupied Sites Applied to the Solid Electrolyte Cubic Li7La3Zr2O12.

To date, experimental and theoretical works have been unable to uncover the ground-state configuration of the solid electrolyte cubic Li7La3Zr2O12 (c-LLZO). Computational studies rely on an initial low-energy structure as a reference point. Here, we present a methodology for identifying energetically favorable configurations of c-LLZO for a crystallographically predicted structure. We begin by eliminating structures that involve overlapping Li atoms based on nearest neighbor counts. We further reduce the configuration space by eliminating symmetry images from all remaining structures. Then, we perform a machine learning-based energetic ordering of all remaining structures. By considering the geometrical constraints that emerge from this methodology, we determine that a large portion of previously reported structures may not be feasible or stable. The method developed here could be extended to other ion conductors. We provide a database containing all of the generated structures with the aim of improving accuracy and reproducibility in future c-LLZO research.

Does a Machine-Learned Potential Perform Better Than an Optimally Tuned Traditional Force Field? A Case Study on Fluorohydrins.

We present a comparative study that evaluates the performance of a machine learning potential (ANI-2x), a conventional force field (GAFF), and an optimally tuned GAFF-like force field in the modeling of a set of 10 γ-fluorohydrins that exhibit a complex interplay between intra- and intermolecular interactions in determining conformer stability. To benchmark the performance of each molecular model, we evaluated their energetic, geometric, and sampling accuracies relative to quantum-mechanical data. This benchmark involved conformational analysis both in the gas phase and chloroform solution. We also assessed the performance of the aforementioned molecular models in estimating nuclear spin-spin coupling constants by comparing their predictions to experimental data available in chloroform. The results and discussion presented in this study demonstrate that ANI-2x tends to predict stronger-than-expected hydrogen bonding and overstabilize global minima and shows problems related to inadequate description of dispersion interactions. Furthermore, while ANI-2x is a viable model for modeling in the gas phase, conventional force fields still play an important role, especially for condensed-phase simulations. Overall, this study highlights the strengths and weaknesses of each model, providing guidelines for the use and future development of force fields and machine learning potentials.

Computing Viscosities of Mixtures of Ester-Based Lubricants at Different Temperatures.

Synthetic esters are used as lubricants for applications at high temperatures, but their development can be a trial and error process. In this context, molecular dynamics simulations could be used as a tool to investigate the properties of new lubricants, in particular viscosity. We employ nonequilibrium molecular dynamics (NEMD) simulations to predict bulk Newtonian viscosities of a set of mixtures of two esters, di(2-ethylhexyl) sebacate (DEHS) and di(2-ethylhexyl) adipate (DEHA) at 293 and 343 K as well as equilibrium molecular dynamics (EMD) and NEMD at 393 K and compare these to experimental measurements. The simulations predict mixture densities within 5% of the experimental values, and we are able to retrieve between 99% and 75% of the experimental viscosities for all ranges of temperature. Experimental viscosities show a linear trend which we are able to capture using NEMD at low temperature and EMD at high temperature. Our work shows that, using EMD and NEMD simulations, and the workflows we developed, we can obtain reliable estimates of the viscosities of mixtures of industrially relevant ester-based lubricants at different temperatures.

Atomistic simulations on the carbidisation processes in Pd nanoparticles.

The formation of interstitial PdC x nanoparticles (NPs) is investigated through DFT calculations. Insights on the mechanisms of carbidisation are obtained whilst the material's behaviour under conditions of increasing C-concentration is examined. Incorporation of C atoms in the Pd octahedral interstitial sites is occurring through the [111] facet with an activation energy barrier of 19.3-35.7 kJ mol-1 whilst migration through the [100] facet corresponds to higher activation energy barriers of 124.5-127.4 kJ mol-1. Furthermore, interstitial-type diffusion shows that C will preferentially migrate and reside at the octahedral interstitial sites in the subsurface region with limited mobility towards the core of the NP. For low C-concentrations, migration from the surface into the interstitial sites of the NPs is thermodynamically favored, resulting in the formation of interstitial carbide. Carbidisation reaction energies are exothermic up to 11-14% of C-concentration and slightly vary depending on the shape of the structure. The reaction mechanisms turn to endothermic for higher concentration levels showing that C will preferentially reside on the surface making the interstitial carbide formation unfavorable. As experimentally observed, our simulations confirm that there is a maximum concentration of C in Pd carbide NPs opening the way for further computational investigations on the activity of Pd carbides in directed catalysis.

Using molecular dynamics to simulate realistic structures of nitrocellulose of different nitration levels.

Nitrocellulose is a reactive derivative of cellulose, one of the most commonly occurring natural materials. Nitration of cellulose decreases the stability of the structure, meaning less is understood about its structure and reactions. Although cellulose is often found in fully crystalline forms, nitrocellulose is more commonly paracrystalline, or amorphous. We present a protocol based on molecular dynamics simulations for creating realistic structures of nitrocellulose, particularly focusing on the crystallinity of the systems being created. We will also provide a detailed analysis of the geometric and dynamical parameters used to quantify the degree of crystallinity for the structures created here, with nitration levels varying from 0-14.14 wt% nitrogen content. Paracrystalline cellulose was not created using the protocol designed here, although it was found that the more nitrated a nitrocellulose system, the more the structure tends to paracrystallinity. This is due to a decrease in the number of hydrogen bonds present, and an increase in the size of the functional groups pushing the chains apart and weakening the interactions between the chains of the structure. The structures created are representative of realistic systems, which in the future will be able to be used to build further understanding of long-term storage of nitrocellulose.

BRD4: quantum mechanical protein-ligand binding free energies using the full-protein DFT-based QM-PBSA method.

Fully quantum mechanical approaches to calculating protein-ligand free energies of binding have the potential to reduce empiricism and explicitly account for all physical interactions responsible for protein-ligand binding. In this study, we show a realistic test of the linear-scaling DFT-based QM-PBSA method to estimate quantum mechanical protein-ligand binding free energies for a set of ligands binding to the pharmaceutical drug-target bromodomain containing protein 4 (BRD4). We show that quantum mechanical QM-PBSA is a significant improvement over traditional MM-PBSA in terms of accuracy against experiment and ligand rank ordering and that the quantum and classical binding energies are converged to a similar degree. We test the interaction entropy and normal mode entropy correction terms to QM- and MM-PBSA.

Corrigendum to "Atomistic level characterisation of ssDNA translocation through the E. coli proteins CsgG and CsgF for nanopore sequencing" [Comput. Struct. Biotechnol. J. 19 (2021) 6417-6430].

[This corrects the article DOI: 10.1016/j.csbj.2021.11.014.].

Understanding Adsorption of Organics on Pt(111) in the Aqueous Phase: Insights from DFT Based Implicit Solvent and Statistical Thermodynamics Models.

Adsorption of organics in the aqueous phase is an area which is experimentally difficult to measure, while computational techniques require extensive configurational sampling of the solvent and adsorbate. This is exceedingly computationally demanding, which excludes its routine use. If implicit solvent could be applied instead, this would dramatically reduce the computational cost as configurational sampling of solvent is not needed. Here, using statistical thermodynamic arguments and DFT calculations with implicit solvent models, we show that semiquantitative values for the free energy and entropy change of adsorption in the aqueous phase (ΔGadssolv and ΔSadssolv) for small organics can be calculated, for a range of coverages. We parametrize the soft sphere based solute dielectric cavity to an approximated free energy of solvation for a single Pt atom at the (111) facet, forming upper and lower bounds based on the entropy of water at the aqueous metal interface (ΔGsolv(Pt) = -4.35 to -7.18 kJ mol-1). This captures the decrease in ΔGadssolv compared to the free energy of adsorption in the vacuum phase (ΔGadsvac), while solvent models with electron density based cavities fail to do so. For a range of oxygenated aromatics, the adsorption energetics using horizontal gas phase geometries significantly overestimate ΔGadssolv compared to experiment by ∼100 kJ mol-1, but they agree with ab initio MD simulations using similar geometries. This suggests oxygenated aromatic compounds adsorb perpendicular to the metallic surface, while the ΔGadssolv for vertical geometries of furfural and cyclohexanol agree to within 20 kJ mol-1 of experimental studies. The proposed techniques provide an inexpensive toolset for validation and prediction of adsorption energetics on solvated metallic surfaces, which could be further validated by the future availability of more experimental measurements for the aqueous entropy/free energy of adsorption.

Energy decomposition analysis method for metallic systems.

In this work, we present the first extension of an energy decomposition analysis (EDA) method to metallic systems. We extend the theory of our Hybrid Absolutely Localized Molecular Orbitals (HALMO) EDA to take into account that molecular orbitals in metallic systems are partially occupied, which is done by weighted orthogonalization (WO) of the molecular orbitals using their associated fractional occupancies as weights in the construction of the projection operators. These operators are needed for the self-consistent field for molecular interaction (SCF MI) computation of the polarization-energy contribution to the interaction. The method gives more weight to orbitals that have higher occupancies and treats each fragment as metallic. The resulting HALMO EDA for metallic systems naturally reduces to the insulator version and produces the same results when applied to an insulating system. We present the theory and implementation of our new approach, and we demonstrate it with sample calculations of relevance to industrial materials. This work provides a new EDA paradigm and tool for the study and analysis of interactions in metallic systems within large-scale DFT calculations.

Atomistic level characterisation of ssDNA translocation through the E. coli proteins CsgG and CsgF for nanopore sequencing.

Two proteins of the Escherichia coli membrane protein complex, CsgG and CsgF, are studied as proteinaceous nanopores for DNA sequencing. It is highly desirable to control the DNA as it moves through the pores, this requires characterisation of DNA translocation and subsequent optimization of the pores. In order to inform protein engineering to improve the pores, we have conducted a series of molecular dynamics simulations to characterise the mechanical strength and conformational dynamics of CsgG and the CsgG-CsgF complex and how these impact ssDNA, water and ion movement. We find that the barrel of CsgG is more susceptible to damage from external electric fields compared to the protein vestibule. Furthermore, the presence of CsgF within the CsgG-CsgF complex enables the complex to withstand higher electric fields. We find that the eyelet loops of CsgG play a key role in both slowing the translocation rate of DNA and modulating the conductance of the pore. CsgF also impacts the DNA translocation rate, but to a lesser degree than CsgG.

Massively parallel linear-scaling Hartree-Fock exchange and hybrid exchange-correlation functionals with plane wave basis set accuracy.

We extend our linear-scaling approach for the calculation of Hartree-Fock exchange energy using localized in situ optimized orbitals [Dziedzic et al., J. Chem. Phys. 139, 214103 (2013)] to leverage massive parallelism. Our approach has been implemented in the onetep (Order-N Electronic Total Energy Package) density functional theory framework, which employs a basis of non-orthogonal generalized Wannier functions (NGWFs) to achieve linear scaling with system size while retaining controllable near-complete-basis-set accuracy. For the calculation of Hartree-Fock exchange, we use a resolution-of-identity approach, where an auxiliary basis set of truncated spherical waves is used to fit products of NGWFs. The fact that the electrostatic potential of spherical waves (SWs) is known analytically, combined with the use of a distance-based cutoff for exchange interactions, leads to a calculation cost that scales linearly with the system size. Our new implementation, which we describe in detail, combines distributed memory parallelism (using the message passing interface) with shared memory parallelism (OpenMP threads) to efficiently utilize numbers of central processing unit cores comparable to, or exceeding, the number of atoms in the system. We show how the use of multiple time-memory trade-offs substantially increases performance, enabling our approach to achieve superlinear strong parallel scaling in many cases and excellent, although sublinear, parallel scaling otherwise. We demonstrate that in scenarios with low available memory, which preclude or limit the use of time-memory trade-offs, the performance degradation of our algorithm is graceful. We show that, crucially, linear scaling with system size is maintained in all cases. We demonstrate the practicability of our approach by performing a set of fully converged production calculations with a hybrid functional on large imogolite nanotubes up to over 1400 atoms. We finish with a brief study of how the employed approximations (exchange cutoff and the quality of the SW basis) affect the calculation walltime and the accuracy of the obtained results.

Generation of Quantum Configurational Ensembles Using Approximate Potentials.

Conformational analysis is of paramount importance in drug design: it is crucial to determine pharmacological properties, understand molecular recognition processes, and characterize the conformations of ligands when unbound. Molecular Mechanics (MM) simulation methods, such as Monte Carlo (MC) and molecular dynamics (MD), are usually employed to generate ensembles of structures due to their ability to extensively sample the conformational space of molecules. The accuracy of these MM-based schemes strongly depends on the functional form of the force field (FF) and its parametrization, components that often hinder their performance. High-level methods, such as ab initio MD, provide reliable structural information but are still too computationally expensive to allow for extensive sampling. Therefore, to overcome these limitations, we present a multilevel MC method that is capable of generating quantum configurational ensembles while keeping the computational cost at a minimum. We show that FF reparametrization is an efficient route to generate FFs that reproduce QM results more closely, which, in turn, can be used as low-cost models to achieve the gold standard QM accuracy. We demonstrate that the MC acceptance rate is strongly correlated with various phase space overlap measurements and that it constitutes a robust metric to evaluate the similarity between the MM and QM levels of theory. As a more advanced application, we present a self-parametrizing version of the algorithm, which combines sampling and FF parametrization in one scheme, and apply the methodology to generate the QM/MM distribution of a ligand in aqueous solution.

Electrochemistry from first-principles in the grand canonical ensemble.

Progress in electrochemical technologies, such as automotive batteries, supercapacitors, and fuel cells, depends greatly on developing improved charged interfaces between electrodes and electrolytes. The rational development of such interfaces can benefit from the atomistic understanding of the materials involved by first-principles quantum mechanical simulations with Density Functional Theory (DFT). However, such simulations are typically performed on the electrode surface in the absence of its electrolyte environment and at constant charge. We have developed a new hybrid computational method combining DFT and the Poisson-Boltzmann equation (P-BE) capable of simulating experimental electrochemistry under potential control in the presence of a solvent and an electrolyte. The charged electrode is represented quantum-mechanically via linear-scaling DFT, which can model nanoscale systems with thousands of atoms and is neutralized by a counter electrolyte charge via the solution of a modified P-BE. Our approach works with the total free energy of the combined multiscale system in a grand canonical ensemble of electrons subject to a constant electrochemical potential. It is calibrated with respect to the reduction potential of common reference electrodes, such as the standard hydrogen electrode and the Li metal electrode, which is used as a reference electrode in Li-ion batteries. Our new method can be used to predict electrochemical properties under constant potential, and we demonstrate this in exemplar simulations of the differential capacitance of few-layer graphene electrodes and the charging of a graphene electrode coupled to a Li metal electrode at different voltages.

Protein-ligand free energies of binding from full-protein DFT calculations: convergence and choice of exchange-correlation functional.

The accurate prediction of protein-ligand binding free energies with tractable computational methods has the potential to revolutionize drug discovery. Modeling the protein-ligand interaction at a quantum mechanical level, instead of relying on empirical classical-mechanics methods, is an important step toward this goal. In this study, we explore the QM-PBSA method to calculate the free energies of binding of seven ligands to the T4-lysozyme L99A/M102Q mutant using linear-scaling density functional theory on the whole protein-ligand complex. By leveraging modern high-performance computing we perform over 2900 full-protein (2600 atoms) DFT calculations providing new insights into the convergence, precision and reproducibility of the QM-PBSA method. We find that even at moderate sampling over 50 snapshots, the convergence of QM-PBSA is similar to traditional MM-PBSA and that the DFT-based energy evaluations are very reproducible. We show that in the QM-PBSA framework, the physically-motivated GGA exchange-correlation functional PBE outperforms the more modern, dispersion-including non-local and meta-GGA-nonlocal functionals VV10 and B97M-rV. Different empirical dispersion corrections perform similarly well but the three-body dispersion term, as included in Grimme's D3 dispersion, is significant and improves results slightly. Inclusion of an entropy correction term sampled over less than 25 snapshots is detrimental while an entropy correction sampled over the same 50 or 100 snapshots as the enthalpies improves the accuracy of the QM-PBSA method. As full-protein DFT calculations can now be performed on modest computational resources our study demonstrates that they can be a useful addition to the toolbox of free energy calculations.

Analysis of DNA interactions and GC content with energy decomposition in large-scale quantum mechanical calculations.

GC content is a contributing factor to the stability of nucleic acids due to hydrogen bonding. More hydrogen bonding generally results in greater stability. Empirical evidence, however, has suggested that the GC content of a nucleic acid is a poor predictor of its stability, implying that there are sequence-dependent interactions besides what its GC content indicates. To examine how much such sequence-dependent interactions affect the interaction energies of double-stranded DNA (dsDNA) molecules, dsDNA molecules of different sequences are generated and examined in silico for variabilities in the interaction energies within each group of dsDNA molecules of the same GC content. Since the amount of hydrogen bonding depends on the GC content, holding the GC content fixed when examining the differences in interaction energies allows sequence-dependent interactions to be isolated. The nature of sequence-dependent interactions is then dissected using energy decomposition analysis (EDA). By using EDA, the components of the interactions that depend on the neighboring base pairs help explain some of the variability in the interaction energies of the dsDNA molecules despite having the same GC content. This work provides a new paradigm and tool for the study and analysis of the distributions of interaction components in dsDNA with the same GC content using EDA within large-scale quantum chemistry calculations.

ParaMol: A Package for Automatic Parameterization of Molecular Mechanics Force Fields.

The ensemble of structures generated by molecular mechanics (MM) simulations is determined by the functional form of the force field employed and its parameterization. For a given functional form, the quality of the parameterization is crucial and will determine how accurately we can compute observable properties from simulations. While accurate force field parameterizations are available for biomolecules, such as proteins or DNA, the parameterization of new molecules, such as drug candidates, is particularly challenging as these may involve functional groups and interactions for which accurate parameters may not be available. Here, in an effort to address this problem, we present ParaMol, a Python package that has a special focus on the parameterization of bonded and nonbonded terms of druglike molecules by fitting to ab initio data. We demonstrate the software by deriving bonded terms' parameters of three widely known drug molecules, viz. aspirin, caffeine, and a norfloxacin analogue, for which we show that, within the constraints of the functional form, the methodologies implemented in ParaMol are able to derive near-ideal parameters. Additionally, we illustrate the best practices to follow when employing specific parameterization routes. We also determine the sensitivity of different fitting data sets, such as relaxed dihedral scans and configurational ensembles, to the parameterization procedure, and discuss the features of the various weighting methods available to weight configurations. Owing to ParaMol's capabilities, we propose that this software can be introduced as a routine step in the protocol normally employed to parameterize druglike molecules for MM simulations.

Translocation of flexible and tensioned ssDNA through in silico designed hydrophobic nanopores with two constrictions.

Protein-inspired nanopores with hydrophobic constriction regions have previously been shown to offer some promise for DNA sequencing. Here we explore a series of pores with two hydrophobic constrictions. The impact of nanopore radius, the nature of residues that define the constriction region and the flexibility of the ssDNA is explored. Our results show that aromatic residues slow down DNA translocation, and in the case of short DNA strands, they cause deviations from a linear DNA conformation. When DNA is under tension, translocation is once again slower when aromatic residues are present in the constriction. However, the lack of flexibility in the DNA backbone provides a narrower window of opportunity for the DNA bases to be retained inside the pore via interaction with the aromatic residues, compared to more flexible strands. Consequently, there is more variability in translocation rates for strands under tension. DNA entry into the pores is correlated to pore width, but no such correlation between width and translocation rate is observed.